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Abstract. We investigate a class of continuum models for the motion of a two- 
dimensional biological group under the influence of nonlocal social interactions. 
The dynamics may be uniquely decomposed into incompressible motion and po- 
tential motion. When the motion is purely incompressible, the model possesses 
solutions which have constant population density and sharp boundaries for all 
time. Numerical simulations of these "swarm patches" reveal rotating mill-like 
swarms with circular cores and spiral arms. The sign of the social interaction 
term determines the direction of the rotation, and the interaction length scale 
affects the degree of spiral formation. When the motion is purely potential, the 
social interaction term has the meaning of repulsion or attraction depending on 
its sign. For the repulsive case, the population spreads and the density profile 
is smoothed. With increasing interaction length scale, the motion becomes more 
convective and experiences slower diffusive smoothing. For the attractive case, 
the population self-organizes into regions of high and low density. The character- 
istic length scale of the density pattern is predicted and confirmed by numerical 
simulations. 



1. Introduction 

1.1. Background 

Examples of collective motion abound in nature. Swarming, schooling, flock- 
ing, and herding have been observed amongst zooplankton, locusts, fish, 
birds, wolves, and other organisms; see for discussions of such groups. 

A remarkable aspect of these aggregations is that individuals move together 
in a coordinated fashion even though interactions between them via sight, 
smell, hearing, or other senses are typically limited to much shorter distances 
than the size of the group. Over the past few decades, mathematical scien- 
tists have begun to tackle the problem of describing how this coordinated 
global structure may arise in biological groups. An overview of modeling 
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Fig. 1. A partial schematic taxonomy of swarming models. 



issues pertaining to swarming is given in 0|, and a partial schematic tax- 
onomy of swarming models, is shown in Figure We now briefly discuss 
this classification of models (noting first that we ignore many other possible 
categorizations, such as deterministic versus stochastic). 

A fruitful approach to modeling swarms has been to treat each indi- 
vidual as a discrete particle. These "individual-based" models have been 
employed in quite a few biological and mathematical studies, including [SJ 
|H|[7||5||5||1U|. These works begin with simple rules of motion for each indi- 
vidual, involving some combination of self-propulsion, random movement, 
and interaction with neighboring organisms. The models typically take the 
form of coupled nonlinear difference equations or ordinary differential equa- 
tions. While numerical simulations of these models have indeed revealed 
collective behavior, a principal disadvantage is that for realistic numbers 
of individuals, analytical results for the collective motion are difficult or 
impossible to obtain. It is worth mentioning that some progress has been 
made in obtaining analytical results for stationary groups. In 10 , a discrete 
model is formulated, and a Lyapunov functional argument is used to success- 
fully predict an equilibrium state of equally-spaced organisms. However, to 
our knowledge, analytical (non-statistical) descriptions of non-equilibrium 
states in discrete swarming models are scarce. 

An alternative approach to studying swarms is to focus on continuum 
equations, which describe relevant quantities as scalar or vector fields. These 
models may be constructed a priori, as in jllU12j . or they may be derived 
by coarse- graining a particle model, as in |5I7|- In general, continuum mod- 
els provide a convenient setting in which to study large populations since 
one may apply machinery from the analysis of partial differential equations 
(PDEs). In the context of swarming, the focus has generally been on models 
in which the population density satisfies a convection-diffusion equation of 
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the form 



p t + V ■ (vp) = V • (DVp). 



(1) 



Here, p(x, t) is the population density, v(x, t) is the velocity field, and x 
is the (one-, two- or three- dimensional) spatial coordinate. This equation 
states that the density is conserved while individuals travel with average 
velocity v. The motion may involve diffusion, whose strength is measured 
by D = D(x,v,p). 

Swarming models of the form of Q may be classified as either dynamic 
or kinematic depending on how the velocity field v is specified. Dynamic 
models couple to Q an equation for the velocity field, such as 



This momentum equation is analogous to Newton's second law. The left- 
hand side is the material (or convective) derivative, i.e. the time derivative 
in a reference frame moving with the velocity field v. The right-hand side 
is simply a sum of forces. The force f represents the self-propulsion of in- 
dividuals, and k is a nonlocal force due interactions with other members 
of the population, to be discussed momentarily. The remaining terms on 
the right side of 10) represent a "viscosity" with strength proportional to 
v, and an external (environmental) force F ext . An example of a dynamic 
model for swarming may be found in 7 . In that work, v — 0, F ext = and 
f (v) = av/|v| — /3v, so that in the absence of social interactions, individuals 
experience a self-propulsion of strength a in their direction of motion and a 
frictional drag of strength j3. In contrast, kinematic models, such as the one 
in ^2], describe the motion of bodies without consideration to the forces 
acting upon them. That is to say, the velocity does not satisfy a momentum 
equation, but rather is simply a functional of the population density, i. e. 



The functional V may include effects like those captured in J5J), such as 
self-propulsion, social interactions, and environmental influence. 

The essence of the swarming phenomenon is the presence of social in- 
teractions. For the velocity equation J2J), these interactions are represented 
by the term k(p,v), and for ©, they are contained in the functional V(/?). 
The social interaction terms might describe effects such as attraction or re- 
pulsion between individuals sufficiently close to each other, or the tendency 
of individuals to orient themselves similarly to their neighbors. Within the 
context of continuum models, then, the social term takes the form of an 
integral operator (most often of convolution form) and the governing equa- 
tions are actually partial integrodifferential equations (PIDEs). Continuum 
models of this form have been studied, for instance, in [7lll2| . 

One challenge associated with continuum models has been the difficulty 
of obtaining biologically realistic swarm solutions, namely, solutions with 
sharp boundaries, relatively constant internal population densities, and long 
lifetimes. For swarms in one spatial dimension, some progress has been made 



v t + v • Vv = f (v) - k(p, v) + t/V 2 v + F ext . 



(2) 



v = V(p). 



(3) 
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in |12j . which also contains extensive background and an associated litera- 
ture review on this issue. We believe that a related issue is the dimensionality 
of the model. Most continuum swarm models, such as those in jl lU12Ul.'{j . 
have only been investigated in one spatial dimension. We expect swarming 
dynamics in higher dimensional models to be qualitatively different since 
that case allows for organisms to vary their orientations continuously, as in 
the "mill" or "vortex" states that have been observed in fish schools, ant 
colonies and other groups |I4lll5IIIfilll7lll8| ). 

One of our goals in this paper is to highlight a difference in the swarming 
problem for one and two dimensions. The possibility of rotational motion 
in two spatial dimensions allows for cohesive swarms with infinite lifetimes, 
even when the velocity rule does not include local drift of organisms. An- 
other goal is to demonstrate a natural way of classifying swarming dynam- 
ics in two spatial dimensions, namely by using the Hodge Decomposition 
Theorem. Our final goal is to make connections between the properties of 
the social interaction terms (for instance, their associated length scales and 
signs) and the large-scale dynamics of the population. 

In the following subsection, we mention some results for constant density 
traveling band solutions of a class of one-dimensional swarm models. These 
are presented for contrast with the two-dimensional case, the study of which 
constitutes the bulk of this paper. In Subsection we state our primary 
results and outline the remainder of the paper. 

1.2. Elementary results for one- dimensional swarms 

A detailed investigation of a swarming model in one spatial dimension has 
been carried out in |12) . In that work, the population density is assumed to 
obey Q). The kinematic velocity rule is 



Here, the first term represents a local density-dependent drift. The remain- 
ing terms are nonlocal components describing attraction and repulsion, with 
the asterisk operator having the meaning of convolution. Note that the re- 
pulsive effects are higher order in the the population density than are the 
attractive ones. The interaction kernel K is odd, piecewise constant, and 
has compact support. It is given by 



where d is an interaction length scale parameter which may be freely chosen. 

Analysis and numerical simulations in |I2) reveal that the model sup- 
ports swarm solutions with biologically realistic characteristics, namely a 
nearly constant internal population density and sharp edges. For density 
independent diffusion, the cohesive swarm has an exponentially long life- 
time before the population is lost through "tails" in the density profile. For 



v = a e p + (A a - A r p)(K * p) 



(4) 
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the case of small density-dependent diffusion, the model has true "traveling 
band" solutions which have compact support. In either case, the cohesive 
motion of the swarm is achieved by an effective cancellation of the social 
interactions. The internal density of the swarm is precisely that at which the 
attractive and repulsive effects cancel each other, so that the only remain- 
ing component of the velocity is local drift. We mention this for contrast 
with the results to be presented later in this paper, which demonstrate a 
nonlocal, i.e. cooperative, means by which a constant density swarm may 
move cohesively. 

We now mention some simple existence and uniqueness results for one- 
dimensional swarms with no diffusion. The population density p satisfies 
the convection equation 

Pt + d x {vp) = 0. (6) 

We will show how in one-dimension, realistic velocity rules which are purely 
nonlocal cannot lead to a constant-speed translation of the population, and 
thus cohesive swarms cannot be maintained. Again, these results are pre- 
sented for contrast with the two-dimensional results given later in this paper. 

Since we are interested in making statements about constant density 
swarms with sharp boundaries, we make the constant density traveling band 
(CDTB) ansatz 

p(x,t) = p a W L (x - ct) (7) 
W L (x - ct)v(x,t) = cW L (x - ct). (8) 

Here, po is the constant population density, c is the speed of the traveling 
band, and Wl(x) is the window function defined without loss of generality 
as 

(J -Ml (9) 

I otherwise. 

The ansatz Q - © automatically satisfies the governing equation (JSJ) in 
S7l, the support of p. Note that we have not placed any restrictions on 
the velocity field outside the support of p since the velocity in unpopulated 
areas is irrelevant to the propagation of the swarm. 

We must also consider an equation defining the velocity field. For con- 
trast, we will consider two velocity rules, each of which may be written as 
a (degenerate) version of the generalized kinematic velocity rule 

v(x,t) = F{p) + G^p) [Kx * Hi(j>)] + G 2 {p) [K 2 * H 2 {p)] (10) 

This is a generalization of Q from F is & functional which captures 
the local density dependence of the velocity. It represents drift velocity of 
organisms, irrespective of social forces. K\ and K 2 are interaction kernels, 
and thus the G[K * H] terms represent nonlocal effects which arise from 
repulsion and attraction between organisms. We will assume that Ki and 
K 2 are integrable. In contrast to ^5j> we wn l further assume that these 
interaction kernels are decreasing in their arguments, so that the influence 
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of the population on a given organism's velocity weakens with distance. We 
allow the functionals F, G, H to be nonlinear for generality. We assume that 
F, G%, G2, H\, H2 are smooth, and that i?i(0) = #2(0) = so that velocities 
from social interactions are induced only by nonzero population. 

The swarm density po and the constant band speed c must satisfy a 
consistency condition via the velocity equation HlOjl in order for the ansatz 
- © to be a solution to ©. 

We first consider the case F = 0, G2 = 0, so that becomes 



v = Gi(p)[Ki*Hi(p)]. (11) 

where the argument of p is understood to be z — x — ct. Equivalently, we 
may obtain a rule of this form by choosing F = 0, Gi = G%, Hi = H2, so 
that 

v = Gi(p)[(Ki+K 2 )*Hi(p)\ (12) 

and both attractive and repulsive effects are represented. Note that for 
these velocity rules, in the absence of interactions, there is no underlying 
(i.e. local) drift velocity. 

Combining JJJ) - (JHJ and (|1 1|) . we obtain the consistency condition 

W L c = W L Gi(p W L )[K 1 *H(p W L )]. (13) 

which may be rewritten as 

z 

c = G°Hi J Ki(C) d( for z e Q L (14) 



where 
and 



G°i = Gi( Po ) (15) 

H°i=Hi{p ) (16) 

By differentiating (|14|l with respect to z and applying the first funda- 
mental theorem of calculus, we see that 

Ki{z) =Ki(z-L) forze^L. (17) 

Thus, Ki satisfying (|17|l must be L-periodic on [— L, L] in order to admit a 
CDTB solution (the structure of Ki outside of [— L, L] is not relevant). We 
will call the set of such kernels T. The logical implication goes the reverse 
way as well, as can be seen from writing down a Fourier series for K\ G Y, 
so that (|14|) is satisfied if and only if K\ e Y. In this case, (|14|l becomes 

c = G°iH°K? (18) 

where 

L 

[ K(()dC (19) 
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There may be families of traveling band solutions parameterized by c which 
bifurcate depending on the structure of the nonlinear functions G\ and Hi . 

However, it is important to realize that the choice K\ S T is not bio- 
logically meaningful, and contradicts our earlier assumption regarding the 
spatial decay of interaction kernels. As discussed above, biologically realistic 
kernels are expected to satisfy dK\jd\z\ < 0, so that for a given individual, 
the effect of other individuals does not increase with distance. However, 
K\ € T cannot satisfy dKi/d\z\ < 0, so at best the kernel would be a con- 
stant, but even this choice is not expected to be a good biological model, 
except perhaps for very small L. 

In contrast to the results just mentioned, we may now consider the case 
G2 = 0, so that HlOfl becomes 

v = F(p) + G 1 {p)[K 1 *H 1 (p)}. (20) 

Equivalently, we may choose K\ = K2, Hi = H^. The velocity rule Q 
in takes this form. Note that now there is a local drift, a self-induced 
contribution to the velocity, which is captured by F. Combining (0) - JEJ 
and 12()p . we obtain the consistency condition 

W L c = W L {F(p W L ) + Gi(poWl) [Ki * H( Po W L )}} (21) 
which may be re-written as 



c = F? + GlH® J K x (C) d( for z E fl L . 

Here, 



(22) 



z-L 



fi=Fi{po) (23) 

and Hi and Gi are given by (|16(l and H15(l . 

There are two cases to consider. If iij'Gj 7^ 0, then ill'L'l) becomes 



F? 



G\H°i 



Z 

J Ki(C)d( forzeX2 L . (24) 

i-L 



This is similar to the previous case. The condition i|24J) may be met only 
for Ki G T, in which case existence and uniqueness of solutions depends 
on the structure of (c - F^)/{G\Hf). For K x £ T, CDTB solutions are not 
possible. 

On the other hand, if iJ"G" = 0, then CDTB are possible for any choice 
of kernel Ki . In this case, the number of possible CDTB solutions depends 
on the number of roots of if"G" for positive p, of which there are expected to 
be a finite number. Looking at the problem from the forward (rather than 
inverse) perspective, for biologically realistic choices of Ki, the velocity 
rule (|20|l will lead to a finite number of CDTB solutions. The densities 
correspond to the solutions pg > of HiG\ = 0, and the wave speeds are 
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given by Fi(pq). Thus, the combination of local and nonlocal velocity terms 
selects particular densities and band speeds, rather than allowing entire 
families of solutions, as in the purely nonlocal case. The allowed CDTB 
densities are those at which the nonlocal interactions disappear. Further, 
since we imagine the total population to be fixed in number, this velocity 
rule actually dictates preferred swarm sizes L. These conclusions are similar 
to those reached in for the particular choice of Fi, G\, H\, K\ given 
by 0. 

1.3. Outline and summary 

The remainder of this paper is devoted to an examination of a nonlocal 
kinematic swarming model in two spatial dimensions. In Section |2| we for- 
mulate an abstract model of animal motion based on simple assumptions. 
We also discuss how the Hodge Decomposition Theorem provides a useful 
way of understanding the two-dimensional motion of the group, namely by 
decomposing it into incompressible motion and potential motion. 

In Section|3we focus on the case of incompressible motion. Since we wish 
to study the motion of a biologically-realistic swarm, we assume the initial 
condition to be a finite group with constant internal population density and 
sharp edges. We show that such a swarm retains these characteristics for all 
time. Numerical simulations demonstrate that the dynamics of the swarm 
are rotational, and that the asymptotic states are vortex-like structures with 
circular cores and a potentially complex arrangement of spiral arms. The 
sign of the social interaction term determines the direction of rotation of the 
swarm, and the characteristic length scale of the interactions determines the 
degree of spiral formation. The spiral states we observed are qualitatively 
similar to the mill states observed in [14I15I16I17I18) ). 

In Section 0] we consider the complementary case of potential motion. 
The sign of the social interaction term determines whether the interac- 
tion represents nonlocal repulsion or attraction. These effects lead, respec- 
tively, to dispersion or aggregation of the population. For the dispersive 
case, shorter interaction length scales result in smoothing of the population 
density profile, while larger interaction length scales lead to motion which 
is more convective. For the case of aggregation, a simple linear stability 
analysis enables us to identify a most unstable wavelength, and thus make 
a prediction about the characteristic length scale of the clumped popula- 
tion distribution that will form. We demonstrate these results by means of 
numerical simulations. 

Finally, we conclude in Section with a brief summary and a discussion 
of directions for future investigation. 

2. A kinematic two-dimensional swarm model 

For the remainder of this paper, we study the dynamics of a two-dimensional 
swarming model. In choosing our model, we make the following assumptions: 
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1. The population density is conserved; birth, death, immigration, and em- 
igration of organisms are negligible on the time scale of the swarming 
dynamics. 

2. The motion of organisms is due solely to social interactions, and thus 
velocities depend nonlocally on the population density (no drift). 

3. Interactions between organisms are pairwisc. 

4. The social interactions are a linear functional of the population density. 

5. The social interactions depend only on the distance between organisms, 
and become weaker with increasing distance. 

Implicit in the second assumption is the supposition that random move- 
ment (e.g. due to fluctuations in the organisms' medium or noise in their 
ability to move) is negligible. The third and fourth, and fifth assumptions 
are made for tractability of the model. The third assumption is made so 
that interaction effects on a given organism will be summable, and this will 
lead to a convolution in our model, similar to the model in |12) . 

In the spirit of the work in |12| . we construct an abstract model, and 
thus we do not incorporate many biological specifics. Our model might be 
interpreted as a one for "flat" (two-dimensional) groups in the absence of 
disturbances such as predators or food sources. Even with the simple as- 
sumptions we have made, the dynamics are complex. As discussed in Section 
03 relaxing some of our assumptions to obtain a more biologically realistic 
model will be an element of our future work. 

Under the assumptions described above, the model takes the form 

Pt + V • (vp) = (25) 
v = y K(|x-y|)p(y)dy = K*p. (26) 

Here, x = (x, y) is the two-dimensional spatial coordinate. Note that l|25[) is 
simply (JTJ with D = 0, and 126(1 is a two-dimensional analog of a degenerate 
case of the velocity rule K is our two-dimensional social interaction 

kernel, which is spatially-decaying and isotropic. 

Since our model includes no drift term, velocities decay in the far field 
and we may apply the Hodge Decomposition Theorem (see, for instance, |19j 1 
This theorem states that a vector field in the plane may be uniquely decom- 
posed into a divergence- free component and a gradient component. That is 
to say, the velocity may be written as 

v = V + V<P, V-^ = 0. (27) 

For smooth vector fields decaying at infinity, the divergence free part has a 
scalar stream function <P satisfying V W = 0. Thus, we can write 

v = V ± !f r + V^>. (28) 

Using an analogy to fluid flow, we may think of \P as a stream function 
for the incompressible part of the flow and as a pressure due to interac- 
tions. For functions with integrable gradients, convolution commutes with 
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derivatives, i.e. (VP) * p — V(P * p), so that for the model 1)25(1 - ((26(1 we 
can directly apply the Hodge decomposition to the interaction kernel K: 

K = V X iV + VP (29) 

where P models the interaction pressure (motion towards and away from 
concentrations of density) and N models additional motion which, as we 
will see, allows for rotation and a cohesive swarm. 

To better understand the model ((25(1 - 1(261) . we separate the dynamics 
into the two cases which we have just discussed, namely incompressible 
motion and potential motion. In the following two sections we study each 
case in turn, and demonstrate how the macroscopic motion of the population 
is affected by the interaction kernel K. 

3. Incompressible motion 

In this case, 

K = V^iV (30) 

so that V • v = 0. Note that this rule has no meaning in one spatial dimen- 
sion, since there is no notion of perpendicular movement. In two dimensions, 
we will see that this type of interaction allows for cohesive movement of the 
swarm without a local drift term. 

The governing equations 1(25(1 - ((26(1 may be written compactly as 

Pt + V • [p(V x iV * p)] = 0. (31) 

We take the scalar interaction function N to be a Gaussian of width d, i. e. 

G d (|x|) = ^ Hx|2/<i2 (32) 

One might include an additional constant prefactor, but this would repre- 
sent a velocity scale and may be removed by rescaling the time variable 
in the equations. The length scale d could also be removed by rescaling, 
in which case the only parameter remaining in the problem would be the 
initial condition. We choose to retain the length scale parameter d since it 
has a clear biological interpretation, and since it is more convenient to vary 
than the length scale of the initial condition. 

A Gaussian interaction was also considered for a linear stability analysis 
in |12j . Other works have used power functions or decaying exponentials \7\ 
I1U| . Our interaction function has somewhat different meaning than the ones 
used in these previous works since it will be applied in two dimensions and 
thus has a rotationally-symmetric structure. We choose Gaussian interac- 
tion functions since they are biologically realistic (in terms of being spatially 
decaying) and because they have convenient mathematical properties such 
as bounded norms and infinite differentiability. We do not intend for this 
choice to be taken too literally as approximating pairwise interactions. Cer- 
tainly other choices of kernels would be equally plausible. Nonetheless, many 
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of our qualitative results will hold true for other classes of smooth, spatially 
decaying interaction functions with normalized integral. 

We will begin by making some general statements about the effect of 
varying the interaction length scale d. For very small values of d, the in- 
teraction function iV resembles a (5-function of strength n. For the limiting 
case d — ► 0, (|31|1 may be written as 



since commutes with convolution and the (5-function acts as the identity 
under convolution. A little algebra reveals that is actually p t — 0, and 
thus the swarm will be stationary. This makes intuitive sense. In the case 
that motion is perpendicular to population gradients in a completely local 
sense, the population density profile cannot change, by construction. Of 
course, in a Lagrangian frame (tracking the coordinates of an individual 
organisms) motion is possible, as long as it is perpendicular to the gradient. 

On the other hand, for very large values of d, N is nearly a constant, 
namely zero. In the formal limit d — > oo, V^N = 0, and once again l|31|l 
becomes p t = 0. This result also makes intuitive sense. In this case, organ- 
isms can sense population gradients infinitely far away, but these gradients 
have no influence on velocity since social interactions are infinitely weak, 
and thus the organisms are stationary. 

For simplicity, and for analogy with the results mentioned in Section fl.2l 
we now focus on constant density solutions of compact support. That is to 
say, we assume that the initial condition is a swarm patch with finite area 
and constant population density pq. By making such a choice, we are not 
modeling the initial formation of a constant-density swarm. Rather, this 
model should be interpreted as a macroscopic description of a swarm in 
which attractive and repulsive forces have already come into balance. We 
will study the subsequent movement of such a swarm. 

We use Green's formula to rewrite (|25|) - (|26|l as an integral over the 
boundary: 



where J? is the support of p, and the boundary dfl is parameterized in a 
clockwise orientation. Here s is the arc length and t is the unit tangent vector 
to the boundary. Following the strategy from fluid dynamics, we adopt 
a Lagrangian framework and track points on the boundary of the swarm 
patch. That is to say, we write down a Langrangian formulation of H34[l 
which will be useful for numerical simulations. Taking a to parameterize 
the boundary of the swarm, we have the equation for z(a,i), the patch 
boundary: 



pt + ttV • \pV x p)\ = 



(33) 




(34) 



dz 




I N(z(a,t) -z(a',t))z a (a',t)da' 



(35) 



where the subscript a indicates a derivative along the boundary. 
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This equation describes a self-deforming curve. From a computational 
standpoint, this formulation is convenient because the dimension of the 
problem has been reduced by one. More importantly, we see that since 
the boundary is a self-deforming curve, the swarm patch retains constant 
internal density and compact support for all time. The philosophy here is 
similar to the contour-dynamics formulation of the two dimensional Euler 
equations |2()j which describe how a fluid region of constant vorticity, or 
vortex patch, evolves in time. The difference is that for the swarm patch 
case, the interaction function N is expected to be spatially decaying in 
order to be biologically meaningful (c/. our modeling assumptions at the 
start of Sectional and our choice in eauationl32|) while for the vortex patch 
problem, N = 4- log l x l- 

Before presenting numerical results of this model, we remark on the 
smoothness of the swarm boundary. For the case of fluid dynamical vortex 
patches mentioned above, solutions of l|35fl with smooth initial data are 
known to stay smooth for all time [2122 . This is the case for our present 
swarm patch problem as well. See the appendix for the sketch of a proof. 

Equation l|35(l may be solved numerically to find the evolution of the 
swarm patch boundary. We now briefly describe our simple numerical al- 
gorithm. An initial swarm patch shape is selected, and the boundary is 
discretized into n nodes. Depending on our initial condition, we take the 
initial number of Lagrangian nodes to be between n = 40 and n = 60. The 
shape of the patch is evolved by using the discretized version of l|35[) . The 
position of each node may be updated by computing its velocity and then 
using a time-stepping rule. We perform the spatial integral in (|35|) using 
Simpson's rule, which is an 0{n 2 ) operation. As a start-up procedure, we 
take three time steps using a fourth order Runge-Kutta method. However, 
since the Runge-Kutta method involves many evaluations of the right-hand 
side of 1|35|) . it is computationally expensive. Thus, we use a fourth-order 
multi-step Adams-Bashforth rule for the remainder of the time steps. We 
take a time step of size At = 0.02. Checks are performed with smaller time 
steps and varying initial discretizations of the swarm patch boundary to 
verify that our solutions are sufficiently well-converged. 

Despite the fact that the boundary stays smooth, numerical simulations 
reveal that it develops complex structure (as we show below) which is also 
a feature of vortex patches [23]. As the swarm patch evolves, it may be 
necessary to re-discretize the boundary in order to have an accurate solu- 
tion (i.e. to have a fine enough mesh to capture new spatial complexity). 
We do so at every time step. Nodes which are adjacent with respect to the 
Lagrangian parameter a are checked for spacing. If the Euclidean distance 
becomes too large, a node is inserted between them using linear interpola- 
tion. Similarly, if nodes become too close together, they are replaced with 
a node whose position is the spatial average of the original ones. While 
this latter step discards detail below a certain length scale, we perform it 
nonetheless so that the total number of nodes does not grow so quickly as to 
make the computation prohibitively slow. Finally, we periodically perform a 
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check to verify that the swarm-patch boundary is not self-intersecting (self- 
intersection of the boundary would break the uniqueness of particle paths 
which the problem must obey). If the boundary is found to self-intersect, 
the simulation is aborted, and must be repeated with a finer threshold of 
spatial detail. 

We note that by symmetry arguments, a rotating disk is an exact solu- 
tion to l|31|l (though the rotation is not solid-body rotation). This is true for 
fluid vortex patches as well; see ^1 f° r a discussion. We assume that N is 
the Gaussian interaction function given by 1)32(1 and calculate the resulting 
velocity of points on the swarm boundary, a circle of radius R. After some 
algebra, we find the velocity |v(i?)|, from which we may compute the period 
of rotation 

T{D) = RrJ\ = fo^ 2 ^ W R *' d2 )Y l ( 36 ) 

where I\ is the modified Bessel function of the first kind of order one. Fig- 
ure|2shows the period of rotation of the boundary, T(d), for a swarm patch 
of radius R = 1 and population density po = I. The exact expression <|36H 
is plotted as a line, while data from a numerical simulation, obtained by 
tracking one of the Lagrangian nodes, is plotted as dots. These results not 
only serve as a check on our algorithm, but demonstrate our previous con- 
clusions about the the behavior of (|31fl in the limits d — > and d — > oo for 
the particular case of a circular patch. 

This example also touches on a connection between the interaction func- 
tion N and the direction of rotation. Since the boundary of the patch has 
a direction associated with its parameterization, the rotation is clockwise 
or counterclockwise according to whether the interaction function N is cho- 
sen to have, respectively, a positive or negative sign. This limitation results 
from our choice of kinematic velocity rule. In the case of a dynamic velocity 
rule, inertial effects would give any initial swarm patch a natural direction 
of rotation. That is to say, for a dynamic rule, the swarm will have the free- 
dom to nucleate a rotational state, for instance, as seen in simulations of 
the model in [7] . In this section, for our kinematic velocity rule, we choose 
N to have a positive sign, and thus swarm patches will always rotate in a 
clockwise manner. 

We now turn to a discussion of the behavior of the model for other 
(noncircular) initial conditions, and for intermediate values of d when some 
nontrivial evolution occurs. We find that for the present case of incompress- 
ible velocity, the dynamics are characterized by an overall rotational motion. 
The solutions at sufficiently long times are vortex-like, as we now illustrate 
with several examples. 

Figure [3| shows the evolution of a swarm patch using the interaction 
function l|32l) with d = 1. The initial boundary of the population is the 
polar curve r(9) = 1 + (1/I0)(cos46>). The square-like initial swarm patch 
experiences clockwise rotation. At time t = 1, the beginnings of spiral arms 
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Fig. 2. Period of rotation T for the boundary of a circular swarm patch under 
the model I3H with the scalar interaction function 1321 . Here, T is shown as a 
function of d, the interaction length scale in 11321 . The line corresponds to the exact 
expression T = 2nR/\v(R)\ given by 13611 . The dots correspond to a numerical 
simulation of the contour dynamics equation 1351 . For this example, the radius of 
the patch is R = 1 and the constant population density is po = 1. 

are visible at the corners of the patch, where the initial curvature was great- 
est. By time t = 3, the spirals have grown longer and the core of the patch is 
becoming circular. This trend continues through the end of the simulation 
at t = 10, at which point the spiral arms have grown even longer and the 
core is nearly a perfect circle. 

The evolution is qualitatively similar even for swarm patches whose ini- 
tial shape is far from circular. For instance, Figure 0] shows the evolution 
of a swarm patch, again with the interaction function (|32l) and d = 1. The 
boundary of the initial shape is an ellipse with a major semiaxis of 1 and 
a minor semiaxis of 0.1. As with the previous example, the swarm patch 
rotates in a clockwise direction. Spiral arms develop at the points furthest 
from the center of the patch, and there is a movement of population density 
towards a developing circular core, which is noticeable at time t — 8 and 
well-defined by t = 10. 

Finally, we comment that a similar evolution occurs even for irregularly 
shaped swarm patches. Figure[5]shows the evolution of such a patch with the 
same interaction function as in the previous two examples. The initial shape 
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Fig. 3. Evolution of a swarm patch under the model ffiip . The boundary of the 
initial shape is a polar curve with radius given by r{0) = 1 + (1/10) (cos 48). The 
constant population density is p — l.The scalar interaction function N is the 
Gaussian given by 132H with interaction length scale d = 1. The constant density 
swarm patch rotates clockwise and develops spiral arms. 
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Time = Time = 6 
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Fig. 4. Evolution of a swarm patch under the model l).'-ilft . The scalar interaction 
function N is the Gaussian given by 1321 with interaction length scale d — 1. The 
constant population density is po = 1. The initial shape is an ellipse with a major 
semiaxis of 1 and a minor semiaxis of 0.1. The patch eventually evolves into a 
circular core with two spiral arms. 
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is generated by the polar function r — f{&) where / is a superposition of 
cosine components with randomly chosen amplitudes and randomly chosen 
low-integer frequencies. As with the previous examples, the patch rotates 
clockwise, developing a circular core and an irregular arrangement of spiral 
arms. 

We have shown in this section that in the case of incompressible mo- 
tion, our simple nonlocal kinematic model has constant density solutions of 
compact support. It was seen directly from 1|31|) that the evolution of any 
initial swarm patch slows for very large or very small values of the inter- 
action length scale d. For intermediate values of d, numerical simulations 
demonstrated that the evolution of these swarm patches is rotational, with 
the direction of motion set a priori by the sign on the interaction func- 
tion N . There is a flux of population towards the rotational center of the 
swarm, where a circular core develops. Spiral arms form at regions of the 
boundary where the curvature is very high. We saw that all of our numerical 
simulations resulted in asymptotic vortex states. 

We close this section by discussing the biological significance of incom- 
pressible velocities. For a constant density swarm in which attractive and 
repulsive effects are in balance, organisms wishing to maintain constant den- 
sity must move with incompressible velocities. Potential velocities, discussed 
in the next section, will lead to variations in population density. Thus, one 
may think of the incompressible velocity terms as those which model the 
aggregate, cooperative dynamics of organisms striving to maintain equal 
spacing. Of course, equal spacing might also be maintained by means of a 
constant local drift, but this is not a cooperative effect, and will not lead to 
the type of vortex-like structure seen here. 



K = VP (37) 
(|26|l may be written compactly as 



4. Potential motion 

In this case, 

so that the model (|25|l - 



p t + V • [p(VP * p)} = 0. 



(38) 



We take the scalar interaction function P — ^fGd where Gd is the Gaussian 
distribution of width d given by (|32|l . In Sectional we made the same choice 
for the scalar interaction function N . In that case, the sign of N simply 
determined the direction of rotation of the swarm. For the present case, we 
will see that the sign of P has a much more dramatic effect on the evolution 
of the population. Specifically, it will determine whether organisms disperse 
or aggregate, as we discuss in the following two subsections. 

4-.1. Dispersion 

Here, we take P = —Gd- Note that (|38() has an analogy to Darcy's law 
for flow in porous media. In fact, in the limit d — > 0, — Gd becomes a 
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Fig. 5. Evolution of a swarm patch under the model l).'-ilft . The scalar interaction 
function N is the Gaussian given by 132H with interaction length scale d = 1. 
The constant population density is po = 1- The initial shape is generated by the 
polar function r — f(9) where / is a superposition of cosine components with 
randomly chosen amplitudes and randomly chosen low-integer frequencies. The 
patch eventually evolves into a circular core with an irregular arrangement of 
spiral arms. We begin the simulation with 84 Lagrangian nodes on the swarm 
patch boundary and end with 485 nodes. 
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(5-function of strength ir and the governing equation (|38Jl is a the porous 
media equation. This is a well-studied PDE which possesses an exact self- 
similar solution, called Barenblatt's solution. For an initial population of 
size Q placed at the origin, Barenblatt's solution is 



A discussion of the porous media equation as it relates to biological disper- 
sal, along with a more general statement of Barenblatt's solution, may be 
found in |24|. In the opposite limit d — > oo, i.e. when social interactions are 
extremely nonlocal, a bit of algebra again reveals that the equation becomes 
simply the steady state pt = 0. The intuitive statement of this limiting case 
is similar to that in the previous section: organisms can sense population 
gradients infinitely far away, but these gradients have no influence on ve- 
locity because the strength of the social interactions is infinitely weak. 

For intermediate values of the interaction length scale d, the population 
density profile experiences diffusion and convection. We may understand 
this better by writing the governing equation l|38|l in an alternate form. 
After some algebra, we obtain 



The first term on the right-hand side of l|40|l is convective, and due to 
the single derivative on G, scales like l/d 4 . In contrast, the second term is 
diffusive, and scales like l/d 6 . Thus, for a given population density profile, 
as d is increased, we expect that convection will be more dominant than 
diffusion. 

We demonstrate the role that the interaction length scale d plays by 
means of numerical simulations. For these examples, we focus on a radially 
symmetric model, so that the density p(r, t) is a function of the radial co- 
ordinate r and the velocity v is as given above. Note that if p is radially 
symmetric and P is also radially symmetric then the velocity field for this 
gradient flow points in the radial direction and is itself radially symmetric. 
We solve the governing equation on the unit disk with boundary condition 
p = on the circumference. We use MacCormack's method, which is second- 
order accurate in space and time; see, for instance, |25j . We use n = 64 grid 
points with timesteps of At = 1 x 10~ 5 — 5 x 1CP 4 . Checks are performed 
with finer meshes in space and time to verify that solutions are sufficiently 
well-converged. 

Two example simulations are shown in FigureEl We choose as a random 
initial condition the function 



Here, f(r) is created by superposing Fourier modes with low integer wave 
numbers and random coefficients. Multiplication by the bracketed combi- 
nation is carried out so that the initial condition decays smoothly towards 




Pt = Vp • (VG d * p) + p{V 2 G d * p) 



(40) 




(41) 
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zero. For the first example, this state is evolved with d = 0.01, so that social 
interactions are only very slightly nonlocal. In this case, (|38|) is nearly the 
porous media equation, and the "bumpy" initial condition quickly smooths 
out and approaches the parabolic profile given by Barenblatt's solution. We 
verify that the numerical solution approaches Barenblatt's solution as fol- 
lows. We fit the numerical solution at time t = 0.006 to Ij39(l . The numerical 
solution is evolved numerically, and the fit Barenblatt's solution is evolved 
analytically. The two are compared again at time t — 0.009. Both curves 
are contained in Figure The curves nearly overlay each other, and the 
maximum error between the two is 3%. 

For contrast, we have taken the same random initial condition and in- 
tegrated it with the more nonlocal interaction length scale d — 0.5. In this 
case, Fourier modes are damped much more slowly, and the bumpy initial 
conditions retains its shape much longer. The motion of the swarm is much 
more convective, and the population is transported away from the origin. 

4-2. Aggregation 

In this case, we take P = Gd- Whereas P was strictly negative in the 
previous subsection, it is now strictly positive, and this change has dramatic 
consequences for the dynamics. Now, the governing equation states that 
velocities are up, rather than down, population gradients (nonlocally) , so 
that the population will tend to form groups. 

We may understand this grouping by means of a linear stability analysis. 
To do so, we consider small perturbations p to a constant density steady 
state pq. Linearizing l|38|) . we obtain 

Pt = -poG d * V 2 /) (42) 

and thus we see that the perturbation obeys a nonlocal backwards heat 
equation. Taking a Fourier ansatz for the perturbation i.e. p = pQe l ^' x+crt ^ , 
we find that the linear growth rate is given by 

a(k) =irp Q k 2 e~ k2 ^ 4d2 \ (43) 

where k = |k|. By computing the critical points of (|43|l we see that the 
most unstable modes are those with wave number k u = 2/d. The growth 
of this most unstable mode provides a mechanism for the clumping of or- 
ganisms. We expect that extremely localized interactions will lead to the 
formation of a larger number of small groups, i.e. a density distribution 
pattern with a small characteristic length scale. On the other hand, more 
nonlocal interactions will result in a smaller number of large groups, i.e. a 
density distribution pattern with a larger characteristic length scale. 

We confirm this prediction by means of numerical simulations. Using a 
pseudospectral Fourier method with 128 modes on each axis, we integrate 
(I38|l on a 2ir x 2ir box with periodic boundary conditions. We choose the 
initial density distribution to be p = 1 plus a small random perturbation 
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Fig. 6. Four snapshots of the time evolution of 138H in radially symmetric ge- 
ometry with a random initial condition and the potential interaction function 
P — —Gd given by 132H . The graphs show the population density p versus the 
radial coordinate r. with (a) d = 0.01 Interactions are very localized, and the 
dynamics are similar to those of the porous media equation. The curve labeled 
"Barenblatt" is a snapshot of the self-similar solution to the porous media equa- 
tion given by 1391 and should be compared to the numerical solution at time 
t — 0.009. See text for details, (b) d — 0.5. For this case of more nonlocal interac- 
tions, the population is convected away from the origin and the smoothing of the 
population density profile occurs more slowly. 
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constructed by superposing low wavenumber (fc < 15) Fourier modes with 
random coefficients. For time-stepping, we initialize with a forward Euler 
step and then use a second-order Adams-Bashforth method. Depending on 
the value of the interaction length scale d, we take time steps of At — 
4 x 10~ 5 — 1 x 10~ 3 . Checks are performed with different numbers of modes 
and different time steps to verify convergence. 

Our results are shown in Figured FiguresEJi shows the initial condition. 
Dark patches correspond to regions of higher density. Figure 0d shows the 
center of the power spectrum of the initial perturbation, which is noisy. 
Figure 0; shows the evolution of the state in Figure [TJi at time t = 0.132 
with interaction length scale d = 0.4. Notice the patches of high population 
density. By the linear stability arguments given above, the characteristic 
wave number of the grouping pattern is predicted to be k u — 2/d = 5. 
Figure 0i shows a blow-up of the center part of the power spectrum of 
the evolution of the perturbation. As predicted, the strongest peaks are 
centered around the circle k = 5. Figures and [7J are analogous to 13; 
and0i, but at time t — 2.74 and with the more nonlocal interaction length 
scale d = 1. In this case, fewer groups form, and they are larger. The most 
unstable wave number from linear analysis is k u — 2, and indeed this is the 
wave number corresponding to the strongest peaks in the power spectrum. 
Finally, we comment that these simulations are not continued for longer 
times because they experience exponential blow up, due to the lack of any 
effects to counterbalance the attractive forces in the model. See the appendix 
for a mathematical discussion of the blow-up. Modifications to the model 
to prevent blow up will be a key aspect of future work, as mentioned in the 
next section. 

5. Conclusions 

Our work on biological groups in two dimensions is in the spirit of the one 
dimensional study in j!2| . The overarching goal in this paper is to make 
specific statements about how social interactions between organisms affect 
the large-scale motion of a biological group. In summary, we formulated 
and studied a simple kinematic continuum model which includes nonlocal, 
spatially-decaying social interactions between individuals. We decomposed 
the dynamics of our model into incompressible motion and potential motion, 
or alternatively, motion perpendicular to population gradients and motion 
along gradients. 

Whereas a constant drift and cancellation of nonlocal effects are neces- 
sary for maintaining a cohesive swarm in one-dimension, this is not the case 
in two dimensions. For the special case of an incompressible kernel, so that 
organisms move perpendicular to population gradients (in a nonlocal sense), 
the equations have constant density solutions of compact support. Through 
numerical simulations beginning with a variety of initial conditions, we 
showed that the dynamics for incompressible interactions are rotational, 
and that swarm patches eventually develop vortex-like structure. This rota- 
tional motion is a cooperative mechanism by which a swarm may maintain 
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Fig. 7. Results from integrating 1381 on a 2n x 27r periodic box with interaction 
function P = Gd given by llM2t . (a) Initial population density, given by p — 1 
plus a small random perturbation, (b) Center of the power spectrum of the initial 
perturbation, (c) Population density at t — 0.132 with interaction length scale 
d = 0.4. Note the formation of small, high-density groups, (d) Center of the power 
spectrum of the perturbation at t — 0.132. The strongest peaks are at k = 5, which 
is the most unstable mode as predicted by linear analysis. Figures (e) and (f) are 
analogous to (c) and (d), but data is taken at t = 2.74 and the longer interaction 
length scale d = 1 is used. In this more nonlocal case, the most unstable wave 
number is k — 2 and the population clumps into fewer, larger groups. 
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motion and cohesion once potential (attractive and repulsive) effects have 
come into balance. The sign of the social interaction term determines the 
direction of rotation, and the social interaction length scale determines the 
degree of macroscopic group movement. The observed asymptotic vortex 
states are intriguingly similar to actual mill vortices seen in biological sys- 
tems P^HTMTTTlFTTIFTHj . 

In contrast, potential kernels model repulsion or attraction between or- 
ganisms. For the repulsive case, the interaction length scale determines a 
balance between diffusive motion and convective motion. Very localized in- 
teractions lead to greater smoothing, while more nonlocal interactions result 
in slower smoothing, but more outward convection motion. For the attrac- 
tive case, the length scale determines a most-unstable mode, the growth of 
which results in the clumping of the population into regions of high and low 
density. 

This work leaves open many possibilities for future research. One route 
would be to conduct fully two-dimensional simulations of biological groups 
under the simultaneous influence of incompressible and potential interac- 
tions. Another would be to relax some of the simplifying modeling assump- 
tions made in Section |3 The focus would be on more complicated velocity 
rules containing both nonlocal and local components, each of which may be 
nonlinear. Finally, future work might consider dynamic, rather than kine- 
matic, velocity rules. These rules would take into account inertial forces 
that might capture "phase changes" in animal group behaviors such as the 
transition from milling to translational motion. Ultimately, a model should 
base interaction rules on specific biological socialization functions of the or- 
ganisms. Unfortunately, field and laboratory data leading to such models is 
very limited. We hope that the general discussion of this paper will help to 
focus further research in this direction. 

A. Appendix 

In this appendix we sketch proofs for two results mentioned in the body 
of this paper, namely regularity of a swarm-patch boundary for the model 
examined in Section |21 and an exponential upper bound for the blow-up of 
the model examined in Section l4.2l 

A. 1. Regularity of the boundary of swarm patches 

In this subsection we discuss the swarm patch model of Section |21 and the 
regularity of the swarm patch boundary. The boundary parametrized by 
z(a,t) is a solution of the integrodifferential equation 

dz 

— =v(z,t), v(x,i)=K*p(x,t), p(x,t)=Xn (44) 
where f2(t) is the interior of the swarm patch. 
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If we consider the problem as an integrodifferential equation for p, 

^+v-Vp = 0, v = K*p (45) 

where K = V^iV for a smooth radial function TV decaying at infinity, then 
the swarm patch is an example of a weak solution of this problem with initial 
condition po G L 1 D L°°(R 2 ). Existence and uniqueness of this problem can 
be proved following the classical theory of Yudovich for vortex patches 
which is also detailed in . Such a discussion is beyond the scope of this 
paper. However, we present some straightforward estimates that can be used 
to prove that the boundary of the patch remains smooth if initially smooth, 
as in the case of vortex patches for which the kernel K is more singular (and 
the proof is correspondingly more difficult). 

If the swarm density persists as the characteristic function of a domain 
f2(t) then it is uniformly bounded in L 1 n L°° for all time. Since v = K * p 
with K G C°° then we have an a priori bound for all derivatives of v, 
D k v < oo for all multi-indices k. Smoothness of the patch boundary then 
follows from the fact that the map z satisfies the ODE l|44|) with initial 
condition 

»|to=*o(a) (46) 

for smooth Zq. Since v is C°° by standard regularity theory for solutions 
of ODEs we see that z itself is smooth. Note that the mapping z can not 
develop a critical point at a later time t because the Lagrangian derivative 
of z, Zq, satisfies the ODE 

^ = Vv | z(Q , t) z Q . (47) 

Since Vv is bounded for all time then z Q remains bounded away from zero 
and infinity if it is so bounded at time zero, by Gronwall's Lemma. 

A. 2. Boundedness of the swarm density for a general velocity rule 

In Scction f4.2l we presented numerical computations that showed p(x, i) can 
exhibit blowup when the convolution kernel K = VP = VGd is positive. 
In the limit as d — > this formally corresponds to a backward time porous 
media equation. 

Here we derive an a priori bound that shows that for a smooth kernel of 
any sign, that the maximum of p is bounded by an exponential in time. Thus 
the blowup seen numerically must be an infinite time blowup, not finite time. 
The bound we derive depends on the L°° norm of the convolution kernel. 
Thus as d — > the bound itself becomes unbounded, as it should because 
we are approaching the ill-posed limit in the positive kernel case. 

In Eulerian coordinates, p satisfies a reaction convection equation 



p t + v • Vp = — pV • v. 



(48) 
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This problem can be transformed to Lagrangian coordinates using the method 
of characteristics. Let X(a, t) denote the solution of the ODE 

^=v(X,t), X| t=0 =a. (49) 
Then in the Lagrangian coordinate a, p satisfies 

^ = -(V-v)| x(0it) p. (50) 
Thus we have a differential inequality for the L°° norm of p, 

|||p|U~<C||V-v||^||p||^. (51) 

Since p is a density, it has an a priori L 1 bound, 

J p(x,t)dx = J p(x,0)dx= ||p||ii. (52) 

Since v = K * p for smooth K, we have 

HV-vlU- < ||V-K||iao|[p[| £ i. (53) 

Combining this with H51|l and the a priori bound on the L 1 norm of p, we 
have 

f t \\p\\ L ~<C\\p\\ L i\\p\\ L ~, C=||V-K|| £ -. (54) 
Gronwall's Lemma then gives 

||p||i- <eajp(C||po|| z it)||polU-. (55) 
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